Here I will explore how DEPOSITED AND SUSPENDED sediment affects the PRESENCE/ABSENCE OF NON-ADVERSE AND ADVERSE EFFECTS of tropical, scleractinean coral of all life-history stages. The database to be used is comprised of data extracted from studies deemed relevant during the systematic literature review process.

The Dataset

Deposited Sediment

Specifically, we will explore the results of 1644 datapoints from 45 studies conducted in 3 oceans on 113 species within 53 genera (shown below only the top 10 species and genera, by number of datapoints).

##      n
## 1 1644
##      Ref   n
## 1   DS01  24
## 2   DS02  16
## 3   DS03  12
## 4   DS04  70
## 5   DS05  10
## 6   DS06   6
## 7   DS07  13
## 8  DS08a 162
## 9  DS08b 113
## 10  DS10  16
## 11  DS11   4
## 12  DS12   7
## 13  DS13 126
## 14  DS15  54
## 15  DS16  46
## 16  DS17  24
## 17  DS18  16
## 18  DS19  62
## 19  DS20   8
## 20  DS21   4
## 21  DS23  51
## 22  DS24   4
## 23  DS25   2
## 24  DS26   6
## 25  DS27  39
## 26  DS28   4
## 27  DS29  40
## 28  DS30  36
## 29  DS31 168
## 30  DS32  16
## 31  DS33   8
## 32  DS34 157
## 33  DS36   2
## 34  DS37   4
## 35  DS38  24
## 36  DS40  20
## 37  DS42   2
## 38  DS43  12
## 39  DS45  24
## 40  DS46  44
## 41  DS48  20
## 42  DS49   6
## 43  DS68  48
## 44  DS69   2
## 45  DS71 112
##        Ref                                    Ref_name
## 1     DS01                   Babcock and Davies (1991)
## 25    DS02                    Babcock and Smith (2000)
## 41    DS03              Bessell-Browne et al. (2017a) 
## 53    DS04                     Duckworth et al. (2017)
## 123   DS05                                Hodel (2007)
## 133   DS06                             Hodgson (1990a)
## 139   DS07                             Hodgson (1990b)
## 152  DS08a               Hodgson (1989) PhD Chapter IV
## 314  DS08b               Hodgson (1989) PhD Chapter IV
## 427   DS10                        Junjie et al. (2014)
## 443   DS11                        Lirman et al. (2008)
## 447   DS12                        Loiola et al. (2013)
## 454   DS13                       Moeller et al. (2017)
## 580   DS15                         Perez et al. (2014)
## 634   DS16                Philipp and Fabricius (2003)
## 680   DS17                               Piniak (2007)
## 704   DS18                     Piniak and Brown (2008)
## 720   DS19                       Ricardo et al. (2017)
## 782   DS20                     Riegl and Branch (1995)
## 790   DS21                              Rogers (1979) 
## 794   DS23                                Selim (2007)
## 845   DS24                      Sheridan et al. (2014)
## 849   DS25                  Shore-Maggio et al. (2018)
## 851   DS26                Sofonia (2006) PhD Chapter 3
## 857   DS27                Sofonia (2006) PhD Chapter 4
## 896   DS28                  Sofonia and Anthony (2008)
## 900   DS29         Stafford-Smith (1990) PhD Chapter 4
## 940   DS30                      Stafford-Smith (1992) 
## 976   DS31            Stafford-Smith and Ormond (1992)
## 1144  DS32                       Stewart et al. (2006)
## 1160  DS33                  Vargas-Angel et al. (2006)
## 1168  DS34                        Weber et al. (2006) 
## 1325  DS36                          Zill et al. (2017)
## 1327  DS37                     Fabricius et al. (2003)
## 1331  DS38                          Goh and Lee (2008)
## 1355  DS40            Abdel-Salam (1989) PhD Chapter 3
## 1375  DS42                         Gowan et al. (2014)
## 1377  DS43                   Peters and Pilson (1985) 
## 1389  DS45                               Rogers (1983)
## 1413  DS46                       Stafford-Smith (1993)
## 1457  DS48  Bessell-Browne et al. (2017) Lab component
## 1477  DS49                             Coffroth (1985)
## 1483  DS68                        Flores et al. (2012)
## 1531  DS69                           Gil et al. (2016)
## 1533  DS71 HDR EOC and CSA Ocean Services, Inc. (2014)
##                                       Ref_name   n
## 1             Abdel-Salam (1989) PhD Chapter 3  20
## 2                    Babcock and Davies (1991)  24
## 3                     Babcock and Smith (2000)  16
## 4   Bessell-Browne et al. (2017) Lab component  20
## 5               Bessell-Browne et al. (2017a)   12
## 6                              Coffroth (1985)   6
## 7                      Duckworth et al. (2017)  70
## 8                      Fabricius et al. (2003)   4
## 9                         Flores et al. (2012)  48
## 10                           Gil et al. (2016)   2
## 11                          Goh and Lee (2008)  24
## 12                         Gowan et al. (2014)   2
## 13 HDR EOC and CSA Ocean Services, Inc. (2014) 112
## 14                                Hodel (2007)  10
## 15               Hodgson (1989) PhD Chapter IV 275
## 16                             Hodgson (1990a)   6
## 17                             Hodgson (1990b)  13
## 18                        Junjie et al. (2014)  16
## 19                        Lirman et al. (2008)   4
## 20                        Loiola et al. (2013)   7
## 21                       Moeller et al. (2017) 126
## 22                         Perez et al. (2014)  54
## 23                   Peters and Pilson (1985)   12
## 24                Philipp and Fabricius (2003)  46
## 25                               Piniak (2007)  24
## 26                     Piniak and Brown (2008)  16
## 27                       Ricardo et al. (2017)  62
## 28                     Riegl and Branch (1995)   8
## 29                              Rogers (1979)    4
## 30                               Rogers (1983)  24
## 31                                Selim (2007)  51
## 32                      Sheridan et al. (2014)   4
## 33                  Shore-Maggio et al. (2018)   2
## 34                Sofonia (2006) PhD Chapter 3   6
## 35                Sofonia (2006) PhD Chapter 4  39
## 36                  Sofonia and Anthony (2008)   4
## 37         Stafford-Smith (1990) PhD Chapter 4  40
## 38                      Stafford-Smith (1992)   36
## 39                       Stafford-Smith (1993)  44
## 40            Stafford-Smith and Ormond (1992) 168
## 41                       Stewart et al. (2006)  16
## 42                  Vargas-Angel et al. (2006)   8
## 43                        Weber et al. (2006)  157
## 44                          Zill et al. (2017)   2
##          Ocean    n
## 1      Pacific 1446
## 2     Atlantic   95
## 3       Indian   63
## 4 Indo-Pacific   40
## Selecting by n
##                           Gsp   n
## 1       Montipora_peltiformis 167
## 2          Acropora_millepora 140
## 3      Pocillopora_damicornis 121
## 4         Leptastrea_purpurea  99
## 5        Porites_lobata/lutea  86
## 6         Acropora_hyacinthus  50
## 7            Leptoria_phrygia  46
## 8          Porites_cylindrica  43
## 9  Montipora_aequituberculata  30
## 10                Porites_rus  30
## Selecting by n
##    Updated_Genus   n
## 1      Montipora 319
## 2       Acropora 258
## 3        Porites 205
## 4    Pocillopora 133
## 5     Leptastrea  99
## 6     Turbinaria  47
## 7       Leptoria  46
## 8         Pavona  37
## 9     Echinopora  31
## 10       Oxypora  31

Suspended Sediment

Specifically, we will explore the results of 661 datapoints from 41 studies conducted in 3 oceans on 29 species within 18 genera (shown below only the top 10 species and genera, by number of datapoints).

##     n
## 1 661
##      Ref   n
## 1   SS06 150
## 2   SS08  48
## 3  SS14c  43
## 4  SS20d  37
## 5  SS16a  32
## 6  SS22c  28
## 7  SS17b  26
## 8  SS17a  24
## 9  SS21a  23
## 10  SS15  20
## 11 SS16c  16
## 12  SS27  16
## 13 SS20c  15
## 14  SS01  12
## 15 SS03a  12
## 16  SS05  12
## 17 SS11a  12
## 18 SS12a  12
## 19 SS22a  10
## 20  SS04   9
## 21 SS12b   8
## 22 SS16b   8
## 23 SS20a   8
## 24 SS22b   8
## 25 SS19a   6
## 26 SS20b   6
## 27 SS21b   6
## 28 SS13a   5
## 29 SS13b   5
## 30 SS13c   5
## 31 SS13d   5
## 32 SS14a   5
## 33 SS19b   4
## 34 SS24a   4
## 35 SS24b   4
## 36  SS07   3
## 37 SS11b   3
## 38 SS11c   3
## 39 SS14b   3
## 40  SS25   3
## 41  SS28   2
##                      Ref_name   n
## 1          Browne et al. 2015 150
## 2         Ricardo et al. 2016  66
## 3         Kendall et al. 1985  56
## 4        Humphrey et al. 2008  51
## 5             Liu et al. 2015  50
## 6          Flores et al. 2012  48
## 7                   Rice 1984  46
## 8         Ricardo et al. 2018  29
## 9        Humanes et al. 2017a  20
## 10       Humanes et al. 2017b  20
## 11         Jokiel et al. 2014  20
## 12               Gilmour 1999  18
## 13        Anthony et al. 2007  16
## 14               Anthony 1999  12
## 15 Anthony and Fabricius 2000  12
## 16         Browne et al. 2014  12
## 17        Ricardo et al. 2015  10
## 18 Bessell-Browne et al. 2017   9
## 19                    Te 1992   8
## 20    Erftemeijer et al. 2012   3
## 21          Te 2001 Chapter 6   3
## 22      Dallmeyer et al. 1982   2
##                           Ocean   n
## 1                       Pacific 392
## 2                        Indian 153
## 3                      Atlantic 104
## 4 Indo-Pacific (e.g. Singapore)  12
## Selecting by n
##                           Gsp   n
## 1          Acropora_millepora 131
## 2             Acropora_tenuis  77
## 3        Acropora_cervicornis  56
## 4           Merulina_ampliata  54
## 5         Pachyseris_speciosa  54
## 6          Platygyra_sinensis  54
## 7           Acropora_muricata  50
## 8  Montipora_aequituberculata  30
## 9         Acropora_digitifera  18
## 10        Acropora_intermedia  16
## Selecting by n
##    Updated_Genus   n
## 1       Acropora 348
## 2       Merulina  54
## 3     Pachyseris  54
## 4      Platygyra  54
## 5      Montipora  46
## 6        Porites  25
## 7    Pocillopora  17
## 8     Goniastrea  12
## 9      Cladocora   8
## 10      Manicina   8
## 11   Solenastrea   8

Exploratory plots

First let’s explore all data from all species for which there is binary data about the presence of ‘non-adverse effects’, ‘adverse effects’, and ‘any mortality’ as a result of exposure to deposited and suspended sediment.

DEFINITIONS:
  • ‘Non-adverse effects’ are defined as hydrostatic inflation, tentacular activity, increased mucus production, and/or congealed sediment.
  • ‘Adverse effects’ are defined as a shift from auto- to hetero-trophy, reduced photosynthetic efficiency, reduced P/R ratio, reduced growth rate, bleaching (or reduced chlorophyll A content), algal overgrowth, smothering, and any degree of tissue/colony mortality. Here we consider only adult corals.
  • ‘Any mortality’ is defined as the presence of small, large, or total necroses of coral’s surface area.

Threshold Analyses with Binary Data

Now let’s calculate thresholds based on the binary data explored above.

DEFINITIONS:
  • ‘LOAEL’, or the ‘lowest observed adverse effect level’ is defined as the lowest dose at which there was an observed adverse effect.
  • ‘NOAEL’, or the ‘no observed adverse effect level’ is defined as the highest dose at which there was NOT an observed adverse effect.
  • ‘Adverse effect’ is defined here as any response of a coral individual, colony, or experimental treatment group that may negatively affect the coral’s fitness and/or survival. These adverse effects may include sublethal physiological changes (e.g., significantly reduced growth or photosynthetic rates when compared to coral in ambient/control conditions), bleaching/tissue loss, and mortality. This definition of adverse effect is independent of its magnitude, so while the affect may have potentially reduced fitness, its magnitude may be sufficiently small that the fitness effect is not measurable.

DEPOSITED SEDIMENT

Adverse effect

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1     1
##   NOAEL
## 1     1

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1   0.5
##   NOAEL
## 1   0.5

Any mortality

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1   4.9
##   NOAEL
## 1   4.4

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1 0.917
##   NOAEL
## 1 0.917

SUSPENDED SEDIMENT

Adverse effect

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  3.22
##   NOAEL
## 1  3.22

NOAEL/LOAEL defined by exposure duration

##        LOAEL
## 1 0.04166667
##       NOAEL
## 1 0.0208333

Any mortality

NOAEL/LOAEL defined by exposure concentration

##   LOAEL
## 1  3.22
##   NOAEL
## 1  3.22

NOAEL/LOAEL defined by exposure duration

##   LOAEL
## 1   0.5
##   NOAEL
## 1   0.5

Logistic Regression Model Fitting

DEPOSITED SEDIMENT

Adverse effect

##     n
## 1 943
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: adverseDS3
## Models:
## modDS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS6:     (1 | Gsp)
## modDS9: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS9:     (1 | Ref)
## modDS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS12:     (1 | Gsp) + (1 | Ref)
## modDS15: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS15:     (1 | Ref_name/Ref)
## modDS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS18:     (1 | Gsp) + (1 | Ref_name)
## modDS21: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS21:     (1 | Gsp) + (1 | Ref_name/Ref)
##         npar    AIC    BIC  logLik deviance    Chisq Df Pr(>Chisq)    
## modDS6     5 1129.5 1153.8 -559.75   1119.5                           
## modDS9     5 1234.8 1259.1 -612.42   1224.8   0.0000  0     1.0000    
## modDS12    6 1078.3 1107.3 -533.13   1066.3 158.5888  1     <2e-16 ***
## modDS15    6 1234.9 1264.0 -611.45   1222.9   0.0000  0     1.0000    
## modDS18    6 1084.3 1113.4 -536.16   1072.3 150.5848  0     <2e-16 ***
## modDS21    7 1080.1 1114.0 -533.04   1066.1   6.2525  1     0.0124 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: adverseDS3
## Models:
## modDS10: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDS10:     Ref)
## modDS11: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDS11:     (1 | Gsp) + (1 | Ref)
## modDS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDS12:     (1 | Gsp) + (1 | Ref)
##         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modDS10    4 1106.5 1125.9 -549.26   1098.5                          
## modDS11    5 1085.9 1110.1 -537.93   1075.9 22.6545  1  1.939e-06 ***
## modDS12    6 1078.3 1107.3 -533.13   1066.3  9.5959  1    0.00195 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: adverseDS3
## Models:
## modDS4: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp)
## modDS7: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Ref)
## modDS10: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDS10:     Ref)
## modDS13: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modDS16: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDS16:     Ref_name)
## modDS19: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDS19:     Ref_name/Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDS4     3 1133.9 1148.4 -563.93   1127.9                         
## modDS7     3 1244.0 1258.6 -619.01   1238.0   0.00  0     1.0000    
## modDS10    4 1106.5 1125.9 -549.26   1098.5 139.50  1     <2e-16 ***
## modDS13    4 1242.3 1261.7 -617.13   1234.3   0.00  0     1.0000    
## modDS16    4 1104.2 1123.6 -548.11   1096.2 138.03  0     <2e-16 ***
## modDS19    5 1106.2 1130.5 -548.11   1096.2   0.00  1     0.9991    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days +  
##     (1 | Gsp) + (1 | Ref)
##    Data: adverseDS3
## 
##      AIC      BIC   logLik deviance df.resid 
##   1085.9   1110.1   -537.9   1075.9      938 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0176 -0.6502 -0.1038  0.5777  3.4105 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 5.151    2.270   
##  Ref    (Intercept) 4.313    2.077   
## Number of obs: 943, groups:  Gsp, 91; Ref, 38
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -1.233170   0.544485  -2.265 0.023523 *  
## Sed_level_stand_mg  0.004003   0.001202   3.330 0.000867 ***
## Sed_exposure_days   0.019949   0.004590   4.346 1.38e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_effect_adverse ~ log10_conc + Sed_exposure_days + (1 |  
##     Gsp) + (1 | Ref)
##    Data: adverseDS3
## 
##      AIC      BIC   logLik deviance df.resid 
##   1067.9   1092.2   -529.0   1057.9      938 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8148 -0.6181 -0.0858  0.5431  3.4637 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 5.470    2.339   
##  Ref    (Intercept) 5.289    2.300   
## Number of obs: 943, groups:  Gsp, 91; Ref, 38
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -3.019371   0.711298  -4.245 2.19e-05 ***
## log10_conc         1.236537   0.239249   5.168 2.36e-07 ***
## Sed_exposure_days  0.021173   0.004644   4.559 5.13e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.804878
##    
## p     0   1
##   0 412 116
##   1  68 347
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.8877
##                    R2m       R2c
## theoretical 0.05186506 0.7779635
## delta       0.04946068 0.7418984

## $Gsp

## 
## $Ref

##  Groups Name        Std.Dev.
##  Gsp    (Intercept) 2.3387  
##  Ref    (Intercept) 2.2998
## [1] 9.971703

##                          Est         LL        UL
## (Intercept)       0.04883195 0.01211253 0.1968671
## log10_conc        3.44366632 2.15460790 5.5039424
## Sed_exposure_days 1.02139871 1.01214419 1.0307379
##                            Est           LL          UL
## (Intercept)       9.563777e-04 3.859093e-05  0.02370138
## log10_conc        1.723998e+01 5.856104e+00 50.75334043
## Sed_exposure_days 1.049961e+00 1.028184e+00  1.07219781

For every 10-fold increase in exposure concentration of deposited sediment, the odds of an adverse effect increase by 3.4 times (95% CI 2.2, 5.5, GLMM z = 5.168, p < 0.0001), while holding exposure duration constant and accounting for variability among studies and species. For every 10-fold increase in exposure duration of deposited sediment, the odds of an adverse effect increase by 1.05 times (95% CI 1.03, 1.07, GLMM z = 4.559, p < 0.0001), while holding exposure concentration constant and accounting for variability among studies and species.

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDS11_log <- predict(modDS11_log, newdata=adverseDS3, type="response")
adverseDS3 <- cbind(adverseDS3, pred_modDS11_log)

adverseDS_plot <- ggplot(data = adverseDS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_effect_adverse)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Adverse effects due to sediment exposure?") +
    scale_x_log10(limits=c(0.01,max(adverseDS3$Sed_level_stand_mg)),
                  breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modDS11_log), inherit.aes=FALSE) +
    theme_bw()

adverseDS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

jvalues <- with(adverseDS3, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    adverseDS3$log10_conc <- j
    predict(modDS11_log, newdata = adverseDS3, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(0,1,2,3)], mean)
## 
# get the means with lower and upper quartiles
plotdat1_2 <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_conc values and convert to data frame
plotdat1_2 <- as.data.frame(cbind(plotdat1_2, jvalues))
plotdat1_2 <- plotdat1_2 %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat1_2) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat1_2)
#approx(plotdat1_2$Sed_conc, plotdat1_2$PredictedProbability, xout=1)
sed <- c(1,5,10,50,100,500,1000)
sed2 <- lapply(sed, function(x) approx(plotdat1_2$Sed_conc, plotdat1_2$PredictedProbability, xout=x))
unlist(sed2)
##            x            y            x            y            x            y 
##    1.0000000    0.1922209    5.0000000    0.3039429   10.0000000    0.3593714 
##            x            y            x            y            x            y 
##   50.0000000    0.4971127  100.0000000    0.5567901  500.0000000    0.6854517 
##            x            y 
## 1000.0000000    0.7335196
# plot average marginal predicted probabilities
ggplot(plotdat1_2, aes(x = log10_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat1_2, aes(x = log10_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat1_2, aes(x = Sed_conc, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
adverseDS_plot2 <- ggplot() +
    geom_point(data = adverseDS3, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_effect_adverse)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of adverse\neffect due to sediment exposure",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.01,max(adverseDS3$Sed_level_stand_mg)), 
                  breaks=c(0.01,0.1,1,10,100,1000,loael),
                  label=c("0.01","0.1","1","10","100","1000",round(loael,digits=1))) +
    geom_ribbon(data = plotdat1_2, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat1_2, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat1_2, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loael, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

adverseDS_plot2

Any mortality

##     n
## 1 719
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: anymortDS2
## Models:
## modDSb6: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb6:     (1 | Gsp)
## modDSb9: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb9:     (1 | Ref)
## modDSb12: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb12:     (1 | Gsp) + (1 | Ref)
## modDSb15: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb15:     (1 | Ref_name/Ref)
## modDSb18: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb18:     (1 | Gsp) + (1 | Ref_name)
## modDSb21: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb21:     (1 | Gsp) + (1 | Ref_name/Ref)
##          npar    AIC    BIC  logLik deviance    Chisq Df Pr(>Chisq)    
## modDSb6     5 758.73 781.62 -374.36   748.73                           
## modDSb9     5 840.41 863.30 -415.20   830.41   0.0000  0   1.000000    
## modDSb12    6 733.43 760.90 -360.72   721.43 108.9770  1  < 2.2e-16 ***
## modDSb15    6 839.84 867.31 -413.92   827.84   0.0000  0   1.000000    
## modDSb18    6 742.10 769.56 -365.05   730.10  97.7451  0  < 2.2e-16 ***
## modDSb21    7 735.05 767.10 -360.53   721.05   9.0432  1   0.002637 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: anymortDS2
## Models:
## modDSb10: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDSb10:     Ref)
## modDSb11: Binary_effect_mortality ~ Sed_level_stand_mg + Sed_exposure_days + 
## modDSb11:     (1 | Gsp) + (1 | Ref)
## modDSb12: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modDSb12:     (1 | Gsp) + (1 | Ref)
##          npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modDSb10    4 768.46 786.77 -380.23   760.46                          
## modDSb11    5 738.10 760.99 -364.05   728.10 32.3569  1  1.283e-08 ***
## modDSb12    6 733.43 760.90 -360.72   721.43  6.6698  1   0.009806 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: anymortDS2
## Models:
## modDSb4: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp)
## modDSb7: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Ref)
## modDSb10: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDSb10:     Ref)
## modDSb13: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Ref_name/Ref)
## modDSb16: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDSb16:     Ref_name)
## modDSb19: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modDSb19:     Ref_name/Ref)
##          npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modDSb4     3 790.99 804.73 -392.50   784.99                         
## modDSb7     3 855.44 869.17 -424.72   849.44  0.000  0     1.0000    
## modDSb10    4 768.46 786.77 -380.23   760.46 88.981  1     <2e-16 ***
## modDSb13    4 852.46 870.77 -422.23   844.46  0.000  0     1.0000    
## modDSb16    4 765.01 783.32 -378.51   757.01 87.445  0     <2e-16 ***
## modDSb19    5 767.01 789.90 -378.51   757.01  0.000  1     0.9997    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_effect_mortality ~ log10_conc + Sed_exposure_days + (1 |  
##     Gsp) + (1 | Ref)
##    Data: anymortDS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    721.3    744.2   -355.6    711.3      714 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3408 -0.4784 -0.0902  0.4327  3.2066 
## 
## Random effects:
##  Groups Name        Variance Std.Dev.
##  Gsp    (Intercept) 7.998    2.828   
##  Ref    (Intercept) 5.484    2.342   
## Number of obs: 719, groups:  Gsp, 81; Ref, 30
## 
## Fixed effects:
##                    Estimate Std. Error z value Pr(>|z|)    
## (Intercept)       -6.749716   1.251060  -5.395 6.84e-08 ***
## log10_conc         2.139224   0.444329   4.814 1.48e-06 ***
## Sed_exposure_days  0.032171   0.006279   5.123 3.00e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.8484006
##    
## p     0   1
##   0 363  56
##   1  53 247
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9317
##                    R2m       R2c
## theoretical 0.09992789 0.8234467
## delta       0.09445315 0.7783326

## $Gsp

## 
## $Ref

##  Groups Name        Std.Dev.
##  Gsp    (Intercept) 2.8280  
##  Ref    (Intercept) 2.3419
## [1] 10.40077

##                           Est           LL          UL
## (Intercept)       0.001171212 0.0001008583  0.01360065
## log10_conc        8.492841038 3.5549391862 20.28961541
## Sed_exposure_days 1.032693755 1.0200618185  1.04548212
##                            Est           LL           UL
## (Intercept)       1.779443e-07 6.283564e-10 5.039204e-05
## log10_conc        1.377919e+02 1.854975e+01 1.023550e+03
## Sed_exposure_days 1.076888e+00 1.046799e+00 1.107843e+00

For every 10-fold increase in exposure concentration of deposited sediment, the odds of any tissue mortality increase by 8.5 times (95% CI 3.6, 20.3; GLMM z = 4.814, p < 0.0001), while holding exposure duration constant and accounting for variability among studies and species. For every 10-fold increase in exposure duration of deposited sediment, the odds of any tissue mortality increase by 1.08 times (95% CI 1.05, 1.11; GLMM z = 5.123, p < 0.0001), while holding exposure concentration constant and accounting for variability among studies and species.

Prediction figure

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modDSb11_log <- predict(modDSb11_log, newdata=anymortDS2, type="response")
anymortDS3 <- cbind(anymortDS2, pred_modDSb11_log)

anymortDS_plot <- ggplot(data = anymortDS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_effect_mortality,
      color = Ref)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Any mortality due to sediment exposure?",
         color = "Study") +
    scale_x_log10(limits=c(0.01,1000),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modDSb11_log), inherit.aes=FALSE) +
    theme_bw()

anymortDS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

jvalues <- with(anymortDS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    anymortDS2$log10_conc <- j
    predict(modDSb11_log, newdata = anymortDS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## 
# get the means with lower and upper quartiles
plotdat2_2 <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_conc values and convert to data frame
plotdat2_2 <- as.data.frame(cbind(plotdat2_2, jvalues))
plotdat2_2 <- plotdat2_2 %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat2_2) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat2_2)
sed <- c(1,5,10,50,100,500,1000)
sed2 <- lapply(sed, function(x) approx(plotdat2_2$Sed_conc, plotdat2_2$PredictedProbability, xout=x))
print(unlist(sed2), digits = 3)
##        x        y        x        y        x        y        x        y 
## 1.00e+00 7.96e-02 5.00e+00 1.92e-01 1.00e+01 2.58e-01 5.00e+01 4.40e-01 
##        x        y        x        y        x        y 
## 1.00e+02 5.21e-01 5.00e+02 6.74e-01 1.00e+03 7.25e-01
# plot average marginal predicted probabilities
ggplot(plotdat2_2, aes(x = log10_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat2_2, aes(x = log10_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat2_2, aes(x = Sed_conc, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
anymortDS_plot2 <- ggplot() +
    geom_point(data = anymortDS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_effect_mortality)) +
    labs(x = expression("Sediment exposure concentration (mg/cm"^"2"*"/day)"), 
         y = "Predicted probability of any tissue\nmortality due to sediment exposure",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.01,1000), breaks=c(0.01,0.1,1,10,100,1000,loael2),
                  label=c("0.01","0.1","1","10","100","1000",round(loael2,digits=1))) +
    geom_ribbon(data = plotdat2_2, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat2_2, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat2_2, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loael2, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

anymortDS_plot2

Suspended Sediment

Adverse Effect

##     n
## 1 423
##      Ref   n
## 1   SS01   8
## 2  SS03a   8
## 3   SS04   6
## 4   SS05   9
## 5   SS06 108
## 6   SS07   2
## 7   SS08  30
## 8  SS11a   8
## 9  SS11b   2
## 10 SS11c   2
## 11 SS12a   9
## 12 SS12b   6
## 13 SS13a   4
## 14 SS13b   4
## 15 SS13c   4
## 16 SS13d   4
## 17 SS14a   4
## 18 SS14b   2
## 19 SS14c  38
## 20 SS17a  20
## 21 SS17b  22
## 22 SS19a   3
## 23 SS19b   2
## 24 SS20a   7
## 25 SS20b   4
## 26 SS20c  13
## 27 SS20d  32
## 28 SS21a  20
## 29 SS21b   3
## 30 SS22a   5
## 31 SS22b   4
## 32 SS22c  14
## 33 SS24a   3
## 34 SS24b   3
## 35  SS25   3
## 36  SS27   6
## 37  SS28   1
##                      Ref_name   n
## 1                Anthony 1999   8
## 2  Anthony and Fabricius 2000   8
## 3         Anthony et al. 2007   6
## 4  Bessell-Browne et al. 2017   6
## 5          Browne et al. 2014   9
## 6          Browne et al. 2015 108
## 7       Dallmeyer et al. 1982   1
## 8     Erftemeijer et al. 2012   2
## 9          Flores et al. 2012  30
## 10               Gilmour 1999  12
## 11       Humanes et al. 2017a  15
## 12       Humanes et al. 2017b  16
## 13       Humphrey et al. 2008  44
## 14            Liu et al. 2015  42
## 15        Ricardo et al. 2015   5
## 16        Ricardo et al. 2016  56
## 17        Ricardo et al. 2018  23
## 18                  Rice 1984  23
## 19                    Te 1992   6
## 20          Te 2001 Chapter 6   3
## Data: adverseSS3
## Models:
## modSS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
## modSS9: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS9:     (1 | Ref)
## modSS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS12:     (1 | Gsp) + (1 | Ref)
## modSS15: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS15:     (1 | Ref_name/Ref)
## modSS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS21:     (1 | Gsp) + (1 | Ref_name/Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modSS6     5 358.90 379.14 -174.45   348.90                         
## modSS9     5 390.53 410.77 -190.27   380.53  0.000  0     1.0000    
## modSS12    6 357.26 381.55 -172.63   345.26 35.270  1  2.871e-09 ***
## modSS15    6 378.86 403.14 -183.43   366.86  0.000  0     1.0000    
## modSS18    6 357.16 381.44 -172.58   345.16 21.701  0  < 2.2e-16 ***
## modSS21    7 357.29 385.62 -171.64   343.29  1.870  1     0.1715    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: adverseSS3
## Models:
## modSS4: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp)
## modSS5: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS5:     (1 | Gsp)
## modSS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
##        npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modSS4    3 393.82 405.97 -193.91   387.82                          
## modSS5    4 362.91 379.10 -177.46   354.91 32.9146  1   9.63e-09 ***
## modSS6    5 358.90 379.14 -174.45   348.90  6.0078  1    0.01424 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: adverseSS3
## Models:
## modSS10: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSS10:     Ref)
## modSS11: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS11:     (1 | Gsp) + (1 | Ref)
## modSS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS12:     (1 | Gsp) + (1 | Ref)
##         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modSS10    4 385.49 401.68 -188.74   377.49                          
## modSS11    5 362.14 382.38 -176.07   352.14 25.3477  1  4.787e-07 ***
## modSS12    6 357.26 381.55 -172.63   345.26  6.8734  1   0.008749 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: adverseSS3
## Models:
## modSS16: Binary_effect_adverse ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSS16:     Ref_name)
## modSS17: Binary_effect_adverse ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS17:     (1 | Gsp) + (1 | Ref_name)
## modSS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
##         npar    AIC    BIC  logLik deviance   Chisq Df Pr(>Chisq)    
## modSS16    4 382.83 399.02 -187.42   374.83                          
## modSS17    5 361.73 381.96 -175.86   351.73 23.1040  1  1.535e-06 ***
## modSS18    6 357.16 381.44 -172.58   345.16  6.5686  1    0.01038 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: adverseSS3
## Models:
## modSS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
## modSS12: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS12:     (1 | Gsp) + (1 | Ref)
## modSS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modSS6     5 358.90 379.14 -174.45   348.90                         
## modSS12    6 357.26 381.55 -172.63   345.26 3.6365  1    0.05652 .  
## modSS18    6 357.16 381.44 -172.58   345.16 0.1061  0    < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: adverseSS3
## Models:
## modSS6: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
## modSS18: Binary_effect_adverse ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## modSS6     5 358.90 379.14 -174.45   348.90                       
## modSS18    6 357.16 381.44 -172.58   345.16 3.7427  1    0.05304 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_effect_adverse ~ log10_conc * Sed_exposure_days + (1 |  
##     Gsp) + (1 | Ref_name)
##    Data: adverseSS3
## 
##      AIC      BIC   logLik deviance df.resid 
##    337.0    361.3   -162.5    325.0      417 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7039 -0.3983 -0.1620 -0.0356  9.8721 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Gsp      (Intercept) 9.4806   3.0791  
##  Ref_name (Intercept) 0.3828   0.6187  
## Number of obs: 423, groups:  Gsp, 26; Ref_name, 20
## 
## Fixed effects:
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -7.698662   1.458821  -5.277 1.31e-07 ***
## log10_conc                    1.996090   0.425708   4.689 2.75e-06 ***
## Sed_exposure_days             0.047097   0.018551   2.539   0.0111 *  
## log10_conc:Sed_exposure_days  0.006267   0.011124   0.563   0.5732    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.8534279
##    
## p     0   1
##   0 304  43
##   1  19  57
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9084
##                   R2m       R2c
## theoretical 0.1372417 0.7842088
## delta       0.1215667 0.6946405

## $Gsp

## 
## $Ref_name

##  Groups   Name        Std.Dev.
##  Gsp      (Intercept) 3.07906 
##  Ref_name (Intercept) 0.61869
## [1] 1.856486

##                                       Est           LL           UL
## (Intercept)                  0.0004534332 2.598605e-05  0.007912004
## log10_conc                   7.3602211982 3.195370e+00 16.953547607
## Sed_exposure_days            1.0482241823 1.010796e+00  1.087038738
## log10_conc:Sed_exposure_days 1.0062864184 9.845841e-01  1.028467094
##                                       Est           LL           UL
## (Intercept)                  2.001417e-08 2.767255e-11 1.447524e-05
## log10_conc                   9.910373e+01 1.451118e+01 6.768261e+02
## Sed_exposure_days            1.114545e+00 1.025033e+00 1.211873e+00
## log10_conc:Sed_exposure_days 1.014534e+00 9.648594e-01 1.066767e+00

For every 10-fold increase in exposure concentration of suspended sediment, the odds of experiencing an adverse effect increase by 7.4 times (95% CI 3.2, 17.0; GLMM z = 4.689, p < 0.0001), while holding constant exposure duration and the interaction between concentration and duration, and while accounting for variability among studies and species. For every 10-fold increase in exposure duration of suspended sediment, the odds of experiencing an adverse effect increase by 1.1 times (95% CI 1.0, 1.2; GLMM z = 2.539, p = 0.011), also while holding constant exposure concentration and the interaction between concentration and duration, and while accounting for variability among studies and species.

Prediction figures

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSS18_log <- predict(modSS18_log, newdata=adverseSS3, type="response")
adverseSS3 <- cbind(adverseSS3, pred_modSS18_log)

adverseSS_plot <- ggplot(data = adverseSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_effect_adverse,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Adverse effect due to sediment exposure",
         color = "Study") +
    scale_x_log10(limits=c(0.1,max(adverseSS3$Sed_level_stand_mg)),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSS18_log), inherit.aes=FALSE) +
    theme_bw()

adverseSS_plot

That plot is confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By sediment exposure concentration
jvalues <- with(adverseSS3, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    adverseSS3$log10_conc <- j
    predict(modSS18_log, newdata = adverseSS3, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat3_2 <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_conc values and convert to data frame
plotdat3_2 <- as.data.frame(cbind(plotdat3_2, jvalues))
plotdat3_2 <- plotdat3_2 %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat3_2) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_conc", "Sed_conc")
#head(plotdat3_2)
sed <- c(1,5,10,50,100,500,1000)
sed2 <- lapply(sed, function(x) approx(plotdat3_2$Sed_conc, plotdat3_2$PredictedProbability, xout=x))
print(round(unlist(sed2), digits=3))
##        x        y        x        y        x        y        x        y 
##    1.000    0.022    5.000    0.068   10.000    0.100   50.000    0.212 
##        x        y        x        y        x        y 
##  100.000    0.280  500.000    0.483 1000.000    0.580
# plot average marginal predicted probabilities
ggplot(plotdat3_2, aes(x = log10_conc)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat3_2, aes(x = log10_conc)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat3_2, aes(x = Sed_conc, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
adverseSS_plot2 <- ggplot() +
    geom_point(data = adverseSS3, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_effect_adverse)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of adverse\neffect due to sediment exposure",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.1,max(adverseSS3$Sed_level_stand_mg)), 
                  breaks=c(0.01,0.1,1,10,100,1000,loaelSS),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelSS,digits=1))) +
    geom_ribbon(data = plotdat3_2, aes(x = Sed_conc, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat3_2, aes(x = Sed_conc, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat3_2, aes(x = Sed_conc, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

adverseSS_plot2

By sediment exposure duration
jvalues <- with(adverseSS3, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    adverseSS3$Sed_exposure_days <- j
    predict(modSS18_log, newdata = adverseSS3, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat3_2b <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat3_2b <- as.data.frame(cbind(plotdat3_2b, jvalues))

# better names and show the first few rows
colnames(plotdat3_2b) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat3_2b)

# plot average marginal predicted probabilities
ggplot(plotdat3_2b, aes(x = Sed_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat3_2b, aes(x = Sed_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
adverseSS_plot3 <- ggplot() +
    geom_point(data = adverseSS3, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_effect_adverse,
      color = Ref)) +
    labs(x = "Sediment exposure duration (days)", 
         y = "Predicted probability of adverse\neffect due to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,max(adverseSS3$Sed_exposure_days))) +
    geom_ribbon(data = plotdat3_2b, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat3_2b, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat3_2b, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

adverseSS_plot3

Any mortality

##     n
## 1 261
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## Data: anymortSS2
## Models:
## modSS6: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS6:     (1 | Gsp)
## modSS9: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS9:     (1 | Ref)
## modSS12: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS12:     (1 | Gsp) + (1 | Ref)
## modSS15: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS15:     (1 | Ref_name/Ref)
## modSS18: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
## modSS21: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS21:     (1 | Gsp) + (1 | Ref_name/Ref)
##         npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## modSS6     5 191.90 209.72 -90.949   181.90                         
## modSS9     5 209.77 227.59 -99.886   199.77  0.000  0          1    
## modSS12    6 193.08 214.47 -90.540   181.08 18.692  1  1.536e-05 ***
## modSS15    6 206.02 227.41 -97.011   194.02  0.000  0          1    
## modSS18    6 186.73 208.12 -87.364   174.73 19.294  0  < 2.2e-16 ***
## modSS21    7 188.73 213.68 -87.364   174.73  0.000  1          1    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Data: anymortSS2
## Models:
## modSS16: Binary_effect_mortality ~ Sed_level_stand_mg + (1 | Gsp) + (1 | 
## modSS16:     Ref_name)
## modSS17: Binary_effect_mortality ~ Sed_level_stand_mg + Sed_exposure_days + 
## modSS17:     (1 | Gsp) + (1 | Ref_name)
## modSS18: Binary_effect_mortality ~ Sed_level_stand_mg * Sed_exposure_days + 
## modSS18:     (1 | Gsp) + (1 | Ref_name)
##         npar    AIC    BIC   logLik deviance  Chisq Df Pr(>Chisq)    
## modSS16    4 216.18 230.44 -104.089   208.18                         
## modSS17    5 195.39 213.21  -92.695   185.39 22.788  1  1.809e-06 ***
## modSS18    6 186.73 208.12  -87.364   174.73 10.661  1   0.001094 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Binary_effect_mortality ~ log10_conc * Sed_exposure_days + (1 |  
##     Gsp) + (1 | Ref_name)
##    Data: anymortSS2
## 
##      AIC      BIC   logLik deviance df.resid 
##    185.1    206.5    -86.6    173.1      255 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6218 -0.3392 -0.1546 -0.0635  4.9901 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Gsp      (Intercept) 4.036    2.009   
##  Ref_name (Intercept) 5.147    2.269   
## Number of obs: 261, groups:  Gsp, 21; Ref_name, 11
## 
## Fixed effects:
##                              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -6.58858    1.87960  -3.505 0.000456 ***
## log10_conc                    0.74576    0.67171   1.110 0.266894    
## Sed_exposure_days             0.02309    0.02624   0.880 0.378980    
## log10_conc:Sed_exposure_days  0.02746    0.01556   1.765 0.077639 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Effect estimates

## [1] 0.8773946
##    
## p     0   1
##   0 212  25
##   1   7  17
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
## Area under the curve: 0.9086
##                   R2m       R2c
## theoretical 0.1683615 0.7806518
## delta       0.1302107 0.6037558

## $Gsp

## 
## $Ref_name

##  Groups   Name        Std.Dev.
##  Gsp      (Intercept) 2.0091  
##  Ref_name (Intercept) 2.2687
## [1] 9.666769

There is no significant relationship between exposure concentration of suspended sediment and the odds of any tissue mortality (GLMM z = 1.110, p = 0.267), nor between exposure duration and the odds of any tissue mortality (GLMM z = 0.880, p = 0.379) after accounting for variability among articles and species.

Prediction figures

#Overlaying glmm results on figure, but prediction line is jagged...
pred_modSSb18_log <- predict(modSSb18_log, newdata=anymortSS2, type="response")
anymortSS3 <- cbind(anymortSS2, pred_modSSb18_log)

anymortSS_plot <- ggplot(data = anymortSS3) +
    geom_point(mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_effect_mortality,
      color = Ref)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Bleaching due to sediment exposure",
         color = "Study") +
    scale_x_log10(limits=c(0.1,max(anymortSS3$Sed_level_stand_mg)),breaks=c(0.01,0.1,1,10,100,1000),
                  label=c("0.01","0.1","1","10","100","1000")) +
    geom_line(aes(x = Sed_level_stand_mg, y = pred_modSSb18_log), inherit.aes=FALSE) +
    theme_bw()

anymortSS_plot

That plot is a bit…confusing to interpret. Let’s see if I can plot the average marginal probability, i.e., the average change in probability of the outcome across the range of the predictor of interest. This is described in some detail at the following, useful website: https://stats.idre.ucla.edu/r/dae/mixed-effects-logistic-regression/

By sediment exposure concentration
jvalues <- with(anymortSS2, seq(from = min(log10_conc), to = max(log10_conc), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    anymortSS2$log10_conc <- j
    predict(modSSb18_log, newdata = anymortSS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat4_2 <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in log10_conc values and convert to data frame
plotdat4_2 <- as.data.frame(cbind(plotdat4_2, jvalues))
plotdat4_2 <- plotdat4_2 %>% mutate(Sed_level_linear=10^jvalues)

# better names and show the first few rows
colnames(plotdat4_2) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "log10_mg_L", "mg_L")
#head(plotdat4_2)
sed <- c(1,5,10,50,100,500,1000)
sed2 <- lapply(sed, function(x) approx(plotdat4_2$mg_L, plotdat4_2$PredictedProbability, xout=x))
print(round(unlist(sed2), digits=3))
##        x        y        x        y        x        y        x        y 
##    1.000    0.019    5.000    0.055   10.000    0.082   50.000    0.170 
##        x        y        x        y        x        y 
##  100.000    0.217  500.000    0.323 1000.000    0.369
# plot average marginal predicted probabilities
ggplot(plotdat4_2, aes(x = log10_mg_L)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat4_2, aes(x = log10_mg_L)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

ggplot(plotdat4_2, aes(x = mg_L, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1)) +
  scale_x_log10()

#Overlaying average marginal predicted probabilities on figure
#by Ref
anymortSS_plot2 <- ggplot() +
    geom_point(data = anymortSS2, mapping = aes(
      x = Sed_level_stand_mg, 
      y = Binary_effect_mortality)) +
    labs(x = "Sediment exposure concentration (mg/L)", 
         y = "Predicted probability of any tissue\nmortality to sediment exposure",
         linetype = "Predicted\nProbability") +
    scale_x_log10(limits=c(0.1,max(anymortSS2$Sed_level_stand_mg)), 
                  breaks=c(0.01,0.1,1,10,100,1000,loaelSS2),
                  label=c("0.01","0.1","1","10","100","1000",round(loaelSS2,digits=1))) +
    geom_ribbon(data = plotdat4_2, aes(x = mg_L, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat4_2, aes(x = mg_L, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat4_2, aes(x = mg_L, y = MedianProbability, 
                                  linetype = "solid")) +
    geom_vline(xintercept=loaelSS2, linetype="dashed", color = "red") +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

anymortSS_plot2

By sediment exposure duration
jvalues <- with(anymortSS2, seq(from = min(Sed_exposure_days), to = max(Sed_exposure_days), length.out = 100))

# calculate predicted probabilities and store in a list
pp <- lapply(jvalues, function(j) {
    anymortSS2$Sed_exposure_days <- j
    predict(modSSb18_log, newdata = anymortSS2, type = "response")
})

# average marginal predicted probability across a few different sediment levels
#sapply(pp[c(1, 20, 40, 60, 80, 100)], mean)
## [1] 0.1797186 0.1960244 0.2142042 0.2333993 0.2535694 0.2746604
# get the means with lower and upper quartiles
plotdat4_2b <- t(sapply(pp, function(x) {
    c(M = mean(x), med = median(x), quantile(x, c(0.25, 0.75)))
}))

# add in Sed_exposure_days values and convert to data frame
plotdat4_2b <- as.data.frame(cbind(plotdat4_2b, jvalues))

# better names and show the first few rows
colnames(plotdat4_2b) <- c("PredictedProbability", "MedianProbability", 
                       "LowerQuantile", "UpperQuantile", "Sed_dur")
#head(plotdat4_2b)

# plot average marginal predicted probabilities
ggplot(plotdat4_2b, aes(x = Sed_dur)) + 
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) +
  ylim(c(0, 1))

ggplot(plotdat4_2b, aes(x = Sed_dur)) + 
  geom_ribbon(aes(ymin = LowerQuantile, ymax = UpperQuantile), alpha = 0.15) +
  geom_line(aes(y = PredictedProbability), size = 2) +
  geom_line(aes(y = MedianProbability), size = 0.5) + 
  ylim(c(0, 1))

#Overlaying average marginal predicted probabilities on figure
#by Ref
anymortSS_plot3 <- ggplot() +
    geom_point(data = anymortSS2, mapping = aes(
      x = Sed_exposure_days, 
      y = Binary_effect_mortality,
      color = Ref)) +
    labs(x = "Sediment exposure duration (days)", 
         y = "Predicted probability of any tissue\nmortality to sediment exposure",
         color = "Binary Data\nby Study",
         linetype = "Predicted\nProbability") +
    scale_x_continuous(limits=c(0,max(anymortSS2$Sed_exposure_days))) +
    geom_ribbon(data = plotdat4_2b, aes(x = Sed_dur, y = PredictedProbability, 
                                    ymin = LowerQuantile, ymax = UpperQuantile), 
                alpha = 0.15) +
    geom_line(data = plotdat4_2b, aes(x = Sed_dur, y = PredictedProbability, 
                                  linetype = "twodash")) +
    geom_line(data = plotdat4_2b, aes(x = Sed_dur, y = MedianProbability, 
                                  linetype = "solid")) +
    theme_bw() +
    scale_linetype_manual(values=c("twodash", "solid"), labels = c("Median","Mean"))

anymortSS_plot3